textbook can be found: https://otexts.org/fpp2/

ts objects

annual freq=1, start=1995
quarterly freq=4, start=c(1995, 2)
monthly freq=12, start=c(1995, 9)
daily freq=(7 or 365.25), start=(1 or c(1995, 234))
weekly freq=52.18, start=c(1995, 23)
hourly freq=(24 or 168 or 8766), start=1
half-hourly freq=(48 or…)

retaildata <- read.csv("https://robjhyndman.com/nyc2018/retail.csv")
mytimeseries <- ts(retaildata[,4], frequency=12, start=c(1982,4))

autoplot(a10)

ggseasonplot(a10)

ggseasonplot(a10, polar=TRUE) + ylab("$ million")

ggsubseriesplot(a10)

beer <- window(ausbeer, start=1992)
autoplot(beer)

ggseasonplot(beer)

ggsubseriesplot(beer)

Lab Session 1

autoplot(mytimeseries)

ggseasonplot(mytimeseries)

ggseasonplot(mytimeseries, polar=TRUE) + ylab("$ million")

ggsubseriesplot(mytimeseries)

Lab Session 2

Seasonality is fixed length, cyclic data have variables length

gglagplot(beer, lags=9)

ggAcf(beer)

elec2 <- window(elec, start=1980)
ggAcf(elec2, lag.max=48)

More correction every 12 months
All lags are positive showing upward trend

autoplot(goog)

ggAcf(goog, lag.max=100)

no seasonality
everything is positive

Session Material

wn <- ts(rnorm(36))
autoplot(wn)

ggAcf(wn)

pigs2 <- window(pigs, start=1990)
autoplot(pigs2)

gglagplot(pigs2)

ggAcf(pigs2)

ggseasonplot(pigs2)

ggsubseriesplot(pigs2)

Lab Session 3

autoplot(bicoal)

gglagplot(bicoal)

ggAcf(bicoal)

autoplot(chicken)

gglagplot(chicken)

ggAcf(chicken)

autoplot(bricksq)

gglagplot(bricksq)

ggAcf(bricksq)

ggseasonplot(bricksq)

ggsubseriesplot(bricksq)

Forecasting Benchmarks

meanf: simply take the mean
naive: use last value
snaive: seasonal naive, take the value of the last same season
rwf: random walk forecast, use drift

meanf(goog, h=20)

##      Point Forecast    Lo 80  Hi 80    Lo 95   Hi 95
## 1001       599.4252 446.0804 752.77 364.7754 834.075
## 1002       599.4252 446.0804 752.77 364.7754 834.075
## 1003       599.4252 446.0804 752.77 364.7754 834.075
## 1004       599.4252 446.0804 752.77 364.7754 834.075
## 1005       599.4252 446.0804 752.77 364.7754 834.075
## 1006       599.4252 446.0804 752.77 364.7754 834.075
## 1007       599.4252 446.0804 752.77 364.7754 834.075
## 1008       599.4252 446.0804 752.77 364.7754 834.075
## 1009       599.4252 446.0804 752.77 364.7754 834.075
## 1010       599.4252 446.0804 752.77 364.7754 834.075
## 1011       599.4252 446.0804 752.77 364.7754 834.075
## 1012       599.4252 446.0804 752.77 364.7754 834.075
## 1013       599.4252 446.0804 752.77 364.7754 834.075
## 1014       599.4252 446.0804 752.77 364.7754 834.075
## 1015       599.4252 446.0804 752.77 364.7754 834.075
## 1016       599.4252 446.0804 752.77 364.7754 834.075
## 1017       599.4252 446.0804 752.77 364.7754 834.075
## 1018       599.4252 446.0804 752.77 364.7754 834.075
## 1019       599.4252 446.0804 752.77 364.7754 834.075
## 1020       599.4252 446.0804 752.77 364.7754 834.075

naive(goog, h=20)

##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 1001         813.67 802.4765 824.8634 796.5511 830.7889
## 1002         813.67 797.8401 829.4999 789.4602 837.8798
## 1003         813.67 794.2824 833.0576 784.0192 843.3208
## 1004         813.67 791.2831 836.0569 779.4322 847.9078
## 1005         813.67 788.6407 838.6993 775.3910 851.9490
## 1006         813.67 786.2517 841.0882 771.7374 855.6026
## 1007         813.67 784.0549 843.2851 768.3776 858.9623
## 1008         813.67 782.0101 845.3298 765.2504 862.0896
## 1009         813.67 780.0896 847.2503 762.3133 865.0267
## 1010         813.67 778.2732 849.0668 759.5353 867.8047
## 1011         813.67 776.5455 850.7945 756.8930 870.4470
## 1012         813.67 774.8947 852.4452 754.3684 872.9716
## 1013         813.67 773.3114 854.0285 751.9469 875.3931
## 1014         813.67 771.7879 855.5520 749.6169 877.7231
## 1015         813.67 770.3179 857.0220 747.3688 879.9712
## 1016         813.67 768.8962 858.4438 745.1944 882.1456
## 1017         813.67 767.5182 859.8218 743.0869 884.2530
## 1018         813.67 766.1802 861.1598 741.0406 886.2993
## 1019         813.67 764.8789 862.4611 739.0504 888.2895
## 1020         813.67 763.6114 863.7286 737.1119 890.2280

snaive(goog, h=20)

##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 1001         813.67 802.4765 824.8634 796.5511 830.7889
## 1002         813.67 797.8401 829.4999 789.4602 837.8798
## 1003         813.67 794.2824 833.0576 784.0192 843.3208
## 1004         813.67 791.2831 836.0569 779.4322 847.9078
## 1005         813.67 788.6407 838.6993 775.3910 851.9490
## 1006         813.67 786.2517 841.0882 771.7374 855.6026
## 1007         813.67 784.0549 843.2851 768.3776 858.9623
## 1008         813.67 782.0101 845.3298 765.2504 862.0896
## 1009         813.67 780.0896 847.2503 762.3133 865.0267
## 1010         813.67 778.2732 849.0668 759.5353 867.8047
## 1011         813.67 776.5455 850.7945 756.8930 870.4470
## 1012         813.67 774.8947 852.4452 754.3684 872.9716
## 1013         813.67 773.3114 854.0285 751.9469 875.3931
## 1014         813.67 771.7879 855.5520 749.6169 877.7231
## 1015         813.67 770.3179 857.0220 747.3688 879.9712
## 1016         813.67 768.8962 858.4438 745.1944 882.1456
## 1017         813.67 767.5182 859.8218 743.0869 884.2530
## 1018         813.67 766.1802 861.1598 741.0406 886.2993
## 1019         813.67 764.8789 862.4611 739.0504 888.2895
## 1020         813.67 763.6114 863.7286 737.1119 890.2280

rwf(goog, drift=TRUE, h=50)

##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 1001       814.0912 802.8996 825.2829 796.9751 831.2073
## 1002       814.5125 798.6773 830.3477 790.2946 838.7304
## 1003       814.9338 795.5300 834.3376 785.2582 844.6093
## 1004       815.3550 792.9383 837.7718 781.0716 849.6385
## 1005       815.7763 790.7011 840.8514 777.4271 854.1254
## 1006       816.1976 788.7154 843.6797 774.1673 858.2278
## 1007       816.6188 786.9200 846.3176 771.1984 862.0393
## 1008       817.0401 785.2749 848.8052 768.4595 865.6207
## 1009       817.4613 783.7526 851.1701 765.9083 869.0144
## 1010       817.8826 782.3329 853.4323 763.5140 872.2512
## 1011       818.3039 781.0005 855.6072 761.2533 875.3544
## 1012       818.7251 779.7438 857.7064 759.1083 878.3419
## 1013       819.1464 778.5533 859.7395 757.0646 881.2281
## 1014       819.5676 777.4214 861.7139 755.1106 884.0247
## 1015       819.9889 776.3419 863.6359 753.2366 886.7412
## 1016       820.4102 775.3095 865.5108 751.4347 889.3856
## 1017       820.8314 774.3199 867.3430 749.6982 891.9647
## 1018       821.2527 773.3692 869.1362 748.0212 894.4842
## 1019       821.6739 772.4542 870.8937 746.3988 896.9491
## 1020       822.0952 771.5720 872.6184 744.8267 899.3638
## 1021       822.5165 770.7202 874.3127 743.3010 901.7320
## 1022       822.9377 769.8966 875.9788 741.8184 904.0571
## 1023       823.3590 769.0993 877.6187 740.3759 906.3420
## 1024       823.7803 768.3265 879.2340 738.9710 908.5895
## 1025       824.2015 767.5766 880.8264 737.6012 910.8019
## 1026       824.6228 766.8483 882.3973 736.2643 912.9812
## 1027       825.0440 766.1403 883.9478 734.9586 915.1295
## 1028       825.4653 765.4515 885.4791 733.6821 917.2485
## 1029       825.8866 764.7808 886.9923 732.4333 919.3398
## 1030       826.3078 764.1272 888.4884 731.2108 921.4048
## 1031       826.7291 763.4900 889.9682 730.0132 923.4450
## 1032       827.1503 762.8682 891.4325 728.8393 925.4614
## 1033       827.5716 762.2611 892.8821 727.6879 927.4553
## 1034       827.9929 761.6682 894.3176 726.5580 929.4278
## 1035       828.4141 761.0886 895.7397 725.4486 931.3797
## 1036       828.8354 760.5218 897.1489 724.3588 933.3119
## 1037       829.2566 759.9674 898.5459 723.2879 935.2254
## 1038       829.6779 759.4247 899.9311 722.2349 937.1209
## 1039       830.0992 758.8933 901.3050 721.1992 938.9991
## 1040       830.5204 758.3728 902.6681 720.1801 940.8608
## 1041       830.9417 757.8626 904.0208 719.1769 942.7065
## 1042       831.3630 757.3625 905.3634 718.1891 944.5368
## 1043       831.7842 756.8721 906.6963 717.2160 946.3524
## 1044       832.2055 756.3910 908.0200 716.2572 948.1537
## 1045       832.6267 755.9188 909.3346 715.3121 949.9413
## 1046       833.0480 755.4554 910.6406 714.3803 951.7157
## 1047       833.4693 755.0003 911.9382 713.4613 953.4772
## 1048       833.8905 754.5533 913.2277 712.5547 955.2263
## 1049       834.3118 754.1142 914.5094 711.6601 956.9634
## 1050       834.7330 753.6826 915.7835 710.7771 958.6890

goog %>% meanf(h=20) %>% autoplot

goog %>% naive(h=20) %>% autoplot

goog %>% snaive(h=20) %>% autoplot

goog %>% rwf(drift=TRUE, h=50) %>% autoplot

Forecasting Residuals

If the fitted values aren’t good, then

goog %>% rwf(drift=TRUE, h=500) %>% fitted() -> z
autoplot(goog, series='Data') + autolayer(z, series='Fitted')

## Warning: Removed 1 rows containing missing values (geom_path).

Residuals \(e_t\) should not be correlated, if they are then information left in the residuals shoudl be used in the forecast

If the mean is not 0 then the forecast is biased

fits <- fitted(naive(goog200))
autoplot(goog200, series='Data') + autolayer(fits, series='Fitted')

## Warning: Removed 1 rows containing missing values (geom_path).

res <- residuals(naive(goog200))
autoplot(res)

gghistogram(res, add.normal=TRUE)

## Warning: Removed 1 rows containing non-finite values (stat_bin).

ggAcf(res)

checkresiduals(naive(goog200))

## 
##  Ljung-Box test
## 
## data:  Residuals from Naive method
## Q* = 11.031, df = 10, p-value = 0.3551
## 
## Model df: 0.   Total lags used: 10

checkresiduals also computes the Ljung-Box test to make sure the residuals are white noise, if p-value is large it is white noise

Lab Session 4

beer <- window(ausbeer, start=1992)
fc <- snaive(beer)
autoplot(fc)

res <- residuals(fc)
autoplot(res)

checkresiduals(fc)

## 
##  Ljung-Box test
## 
## data:  Residuals from Seasonal naive method
## Q* = 32.269, df = 8, p-value = 8.336e-05
## 
## Model df: 0.   Total lags used: 8

This tiny p-value means it’s not white noise and we can do better.

Evaluating Forecast Accuracy

\[ MAPE = 100mean(|e_{T+h}|/y_{T+h}) \]

MAPE doesn’t work well when numbers are close to 0 because of division by \(y_{T+h}\)

MAPE also doesn’t make sense when there is no natura 0, like temperature

Instead use Mean Absolute Scaled Error

\[ MASE = T^{-1} \sum_{t=1}^T |y_t \hat{y_{t|t-1}}| / Q \]

where Q is a stable measure of the scale of the time series

googtrain <- window(goog200, end=180)

googfc1 <- meanf(googtrain, h=20)
googfc2 <- naive(googtrain, h=20)
googfc3 <- rwf(googtrain, drift=TRUE, h=20)

accuracy(googfc1, goog200)

##                         ME     RMSE      MAE        MPE      MAPE
## Training set -1.360278e-14 28.74910 20.06324 -0.4195606  4.582412
## Test set      8.242659e+01 82.89387 82.42659 15.9262898 15.926290
##                   MASE      ACF1 Theil's U
## Training set  5.261155 0.9547241        NA
## Test set     21.614609 0.7822826  20.39543

accuracy(googfc2, goog200)

##                     ME      RMSE       MAE       MPE      MAPE     MASE
## Training set  0.602728  6.404477  3.813467 0.1256394 0.8712766 1.000000
## Test set     16.041900 18.291917 16.041900 3.0762438 3.0762438 4.206645
##                     ACF1 Theil's U
## Training set -0.06213755        NA
## Test set      0.78228260  4.539165

accuracy(googfc3, goog200)

##                         ME      RMSE      MAE         MPE      MAPE
## Training set -2.191158e-14  6.376052 3.887099 -0.01363486 0.8882792
## Test set      9.713257e+00 11.335822 9.713257  1.86155541 1.8615554
##                  MASE        ACF1 Theil's U
## Training set 1.019309 -0.06213755        NA
## Test set     2.547094  0.70229763  2.811493

Lab Session 5

train <- window(mytimeseries, end=c(2010, 12))
test <- window(mytimeseries, start=2011)

autoplot(cbind(Training=train, Test=test))

fcst1 <- snaive(train, h=length(test))

accuracy(fcst1, test)

##                    ME     RMSE      MAE      MPE      MAPE      MASE
## Training set 5.502402 27.74935 19.33784 3.369450 10.447161 1.0000000
## Test set     5.227586 23.40458 18.60000 1.619239  8.172045 0.9618449
##                   ACF1 Theil's U
## Training set 0.8703252        NA
## Test set     0.4544630 0.9518744

checkresiduals(fcst1)

## 
##  Ljung-Box test
## 
## data:  Residuals from Seasonal naive method
## Q* = 1256.8, df = 24, p-value < 2.2e-16
## 
## Model df: 0.   Total lags used: 24

fcst1 %>% autoplot() + autolayer(test)

Time Series Cross-Validation

Rob coined the phrase

rolling training data, all starting at the same point

e <- tsCV(goog200, rwf, drift=TRUE, h=1)
sqrt(mean(e^2, na.rm=TRUE))

## [1] 6.233245

Now with pipes

e <- goog200 %>% tsCV(forecastfunction=rwf, drift=TRUE, h=1)
e^2 %>% mean(na.rm=TRUE) %>% sqrt()

## [1] 6.233245

Now get error for h=1, h=2, h=3, etc

e <- goog200 %>% tsCV(forecastfunction=rwf, drift=TRUE, h=12)
e^2 %>% colMeans(na.rm=TRUE) %>% sqrt()

##       h=1       h=2       h=3       h=4       h=5       h=6       h=7 
##  6.233245  8.634156 10.833826 13.032527 14.962450 16.483554 18.108701 
##       h=8       h=9      h=10      h=11      h=12 
## 19.949726 21.648751 23.237067 24.635515 25.870110

Lab Session 6

e1 <- mytimeseries %>% tsCV(forecastfunction=rwf, drift=TRUE, h=12)
e2 <- mytimeseries %>% tsCV(forecastfunction=meanf, h=12)
e3 <- mytimeseries %>% tsCV(forecastfunction=naive, h=12)
e4 <- mytimeseries %>% tsCV(forecastfunction=snaive, h=12)

e1^2 %>% colMeans(na.rm=TRUE) %>% sqrt()

##      h=1      h=2      h=3      h=4      h=5      h=6      h=7      h=8 
## 20.04797 23.06671 24.18534 25.51110 27.27162 30.09422 29.46610 29.93850 
##      h=9     h=10     h=11     h=12 
## 30.69911 31.68924 31.27128 26.49646

e2^2 %>% colMeans(na.rm=TRUE) %>% sqrt()

##      h=1      h=2      h=3      h=4      h=5      h=6      h=7      h=8 
## 67.49756 67.84488 68.19028 68.53594 68.88044 69.22489 69.56606 69.91056 
##      h=9     h=10     h=11     h=12 
## 70.25265 70.59840 70.94401 71.29212

e3^2 %>% colMeans(na.rm=TRUE) %>% sqrt()

##      h=1      h=2      h=3      h=4      h=5      h=6      h=7      h=8 
## 20.01105 22.98583 24.06746 25.36152 27.08889 29.83772 29.20396 29.64566 
##      h=9     h=10     h=11     h=12 
## 30.36709 31.31045 30.89591 26.30046

e4^2 %>% colMeans(na.rm=TRUE) %>% sqrt()

##      h=1      h=2      h=3      h=4      h=5      h=6      h=7      h=8 
## 26.30046 26.33190 26.36324 26.39421 26.42471 26.45617 26.48818 26.52014 
##      h=9     h=10     h=11     h=12 
## 26.55237 26.58396 26.61484 26.64629

For one-ahead, this gives the overall error for each method

e1 <- mytimeseries %>% tsCV(forecastfunction=rwf, drift=TRUE, h=1)
e2 <- mytimeseries %>% tsCV(forecastfunction=meanf, h=1)
e3 <- mytimeseries %>% tsCV(forecastfunction=naive, h=1)
e4 <- mytimeseries %>% tsCV(forecastfunction=snaive, h=1)

e1^2 %>% mean(na.rm=TRUE) %>% sqrt()

## [1] 20.04797

e2^2 %>% mean(na.rm=TRUE) %>% sqrt()

## [1] 67.49756

e3^2 %>% mean(na.rm=TRUE) %>% sqrt()

## [1] 20.01105

e4^2 %>% mean(na.rm=TRUE) %>% sqrt()

## [1] 26.30046

cbind(RWF=e1^2, Mean=e2^2, Naive=e3^2, SNaive=e4^2) %>% na.omit() %>% colMeans()

##       RWF      Mean     Naive    SNaive 
##  411.4290 4685.8665  410.0438  691.7144

Exponential Smoothing

Simple Exponential Smoothing

Robert Goodall Brown

fc <- ses(oil, h=12)
autoplot(fc)

summary(fc)

## 
## Forecast method: Simple exponential smoothing
## 
## Model Information:
## Simple exponential smoothing 
## 
## Call:
##  ses(y = oil, h = 12) 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
## 
##   Initial states:
##     l = 110.8832 
## 
##   sigma:  49.0507
## 
##      AIC     AICc      BIC 
## 576.1569 576.6903 581.8324 
## 
## Error measures:
##                    ME     RMSE      MAE      MPE     MAPE      MASE
## Training set 8.806147 48.03921 35.52053 2.135668 10.73828 0.9796835
##                   ACF1
## Training set 0.1646214
## 
## Forecasts:
##      Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2014       542.3412 479.4802 605.2022 446.2037 638.4788
## 2015       542.3412 453.4468 631.2356 406.3890 678.2935
## 2016       542.3412 433.4701 651.2124 375.8372 708.8453
## 2017       542.3412 416.6287 668.0537 350.0805 734.6019
## 2018       542.3412 401.7911 682.8914 327.3883 757.2941
## 2019       542.3412 388.3767 696.3057 306.8729 777.8096
## 2020       542.3412 376.0410 708.6415 288.0069 796.6755
## 2021       542.3412 364.5591 720.1233 270.4469 814.2355
## 2022       542.3412 353.7751 730.9074 253.9542 830.7283
## 2023       542.3412 343.5753 741.1072 238.3549 846.3275
## 2024       542.3412 333.8739 750.8085 223.5180 861.1644
## 2025       542.3412 324.6044 760.0781 209.3415 875.3410

summary(fc$model)

## Simple exponential smoothing 
## 
## Call:
##  ses(y = oil, h = 12) 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
## 
##   Initial states:
##     l = 110.8832 
## 
##   sigma:  49.0507
## 
##      AIC     AICc      BIC 
## 576.1569 576.6903 581.8324 
## 
## Training set error measures:
##                    ME     RMSE      MAE      MPE     MAPE      MASE
## Training set 8.806147 48.03921 35.52053 2.135668 10.73828 0.9796835
##                   ACF1
## Training set 0.1646214

Holt’s Linear Trend

Charles Holt (US Navy) (1957)

window(ausair, start=1990, end=2004) %>% 
    holt(h=5, PI=FALSE) %>% 
    summary()

## 
## Forecast method: Holt's method
## 
## Model Information:
## Holt's method 
## 
## Call:
##  holt(y = ., h = 5, PI = FALSE) 
## 
##   Smoothing parameters:
##     alpha = 1e-04 
##     beta  = 1e-04 
## 
##   Initial states:
##     l = 19.1015 
##     b = 1.4307 
## 
##   sigma:  1.9943
## 
##      AIC     AICc      BIC 
## 66.67669 73.34336 70.21694 
## 
## Error measures:
##                       ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -0.09267432 1.707786 1.474578 -0.8251858 5.100699 0.8045133
##                   ACF1
## Training set 0.4475061
## 
## Forecasts:
##      Point Forecast
## 2005       41.99044
## 2006       43.42100
## 2007       44.85155
## 2008       46.28211
## 2009       47.71267

window(ausair, start=1990, end=2004) %>% 
    holt(h=5, PI=FALSE) %>% 
    autoplot()

Damped Trend Method

1980s: Gardner (Houston) & McKinsey (Scotland)

autoplot(livestock)

livestock2 <- livestock %>% window(start=1970, end=2000)
livestock2 %>% holt(h=20, damped=TRUE) %>% autoplot()

Beta=0.3 is pretty big

phi=0.8 is artifically constrained to be the smallest phi can be

Lab Session 7

SES for non-trended, Holt for trended

eggs2 <- window(eggs, start=1900, end=1993)

fc1 <- ses(eggs2, h=100)
fc2 <- holt(eggs2, h=100)
fc3 <- holt(eggs2, h=100, damped=TRUE)

fc1 %>% autoplot()

fc2 %>% autoplot()

fc3 %>% autoplot()

fc1 %>% autoplot(PI=FALSE)

fc2 %>% autoplot(PI=FALSE)

fc3 %>% autoplot(PI=FALSE)

residuals(fc1)^2 %>% mean() %>% sqrt()

## [1] 26.56395

residuals(fc2)^2 %>% mean() %>% sqrt()

## [1] 26.58219

residuals(fc3)^2 %>% mean() %>% sqrt()

## [1] 26.54019

fc1 %>% accuracy()

##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set -2.740622 26.56395 19.38206 -2.872958 10.05648 0.9560879
##                      ACF1
## Training set -0.005983602

fc2 %>% accuracy()

##                      ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.04499087 26.58219 19.18491 -1.142201 9.653791 0.9463626
##                    ACF1
## Training set 0.01348202

fc3 %>% accuracy()

##                     ME     RMSE     MAE       MPE     MAPE      MASE
## Training set -2.891496 26.54019 19.2795 -2.907633 10.01894 0.9510287
##                      ACF1
## Training set -0.003195358

list(SES=fc1, Trend=fc2, DampedTrend=fc3) %>% purrr::map_df(~ accuracy(.x) %>% as.data.frame(), .id='Model')

##         Model          ME     RMSE      MAE       MPE      MAPE      MASE
## 1         SES -2.74062165 26.56395 19.38206 -2.872958 10.056482 0.9560879
## 2       Trend  0.04499087 26.58219 19.18491 -1.142201  9.653791 0.9463626
## 3 DampedTrend -2.89149552 26.54019 19.27950 -2.907633 10.018944 0.9510287
##           ACF1
## 1 -0.005983602
## 2  0.013482023
## 3 -0.003195358

if you care about mean, use RMSE
if you care about median, use MAE or MASE

fc1 %>% checkresiduals()

## 
##  Ljung-Box test
## 
## data:  Residuals from Simple exponential smoothing
## Q* = 10.475, df = 8, p-value = 0.2333
## 
## Model df: 2.   Total lags used: 10

fc2 %>% checkresiduals()

## 
##  Ljung-Box test
## 
## data:  Residuals from Holt's method
## Q* = 11.026, df = 6, p-value = 0.08759
## 
## Model df: 4.   Total lags used: 10

fc3 %>% checkresiduals()

## 
##  Ljung-Box test
## 
## data:  Residuals from Damped Holt's method
## Q* = 10.782, df = 5, p-value = 0.05587
## 
## Model df: 5.   Total lags used: 10

Seasonal Methods

Holt Winters: Holt wrote the pape in 1957, separately Winters (Holt’s students) wrote the code 3 years later

aust <- window(austourists, start=2005)
fc1 <- hw(aust, seasonal='additive')
fc2 <- hw(aust, seasonal='multiplicative')

fc1 %>% autoplot()

fc2 %>% autoplot()

fc1 %>% summary()

## 
## Forecast method: Holt-Winters' additive method
## 
## Model Information:
## Holt-Winters' additive method 
## 
## Call:
##  hw(y = aust, seasonal = "additive") 
## 
##   Smoothing parameters:
##     alpha = 0.3063 
##     beta  = 1e-04 
##     gamma = 0.4263 
## 
##   Initial states:
##     l = 32.2597 
##     b = 0.7014 
##     s=1.3106 -1.6935 -9.3132 9.6962
## 
##   sigma:  1.9494
## 
##      AIC     AICc      BIC 
## 234.4171 239.7112 250.4748 
## 
## Error measures:
##                       ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 0.008115785 1.763305 1.374062 -0.2860248 2.973922 0.4502579
##                     ACF1
## Training set -0.06272507
## 
## Forecasts:
##         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2016 Q1       76.09837 73.60011 78.59664 72.27761 79.91914
## 2016 Q2       51.60333 48.99039 54.21626 47.60718 55.59947
## 2016 Q3       63.96867 61.24582 66.69153 59.80443 68.13292
## 2016 Q4       68.37170 65.54313 71.20027 64.04578 72.69762
## 2017 Q1       78.90404 75.53440 82.27369 73.75061 84.05747
## 2017 Q2       54.40899 50.95325 57.86473 49.12389 59.69409
## 2017 Q3       66.77434 63.23454 70.31414 61.36069 72.18799
## 2017 Q4       71.17737 67.55541 74.79933 65.63806 76.71667

fc2 %>% summary()

## 
## Forecast method: Holt-Winters' multiplicative method
## 
## Model Information:
## Holt-Winters' multiplicative method 
## 
## Call:
##  hw(y = aust, seasonal = "multiplicative") 
## 
##   Smoothing parameters:
##     alpha = 0.4406 
##     beta  = 0.0134 
##     gamma = 0.0023 
## 
##   Initial states:
##     l = 32.4875 
##     b = 0.6974 
##     s=1.0237 0.9618 0.7704 1.2442
## 
##   sigma:  0.0367
## 
##      AIC     AICc      BIC 
## 221.1313 226.4254 237.1890 
## 
## Error measures:
##                      ME     RMSE     MAE           MPE    MAPE      MASE
## Training set 0.09206228 1.575631 1.25496 -0.0006505533 2.70539 0.4112302
##                     ACF1
## Training set -0.07955726
## 
## Forecasts:
##         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## 2016 Q1       80.08894 76.31865 83.85922 74.32278 85.85509
## 2016 Q2       50.15482 47.56655 52.74309 46.19640 54.11324
## 2016 Q3       63.34322 59.80143 66.88502 57.92652 68.75993
## 2016 Q4       68.17810 64.08399 72.27221 61.91670 74.43950
## 2017 Q1       83.80112 78.43079 89.17146 75.58790 92.01434
## 2017 Q2       52.45291 48.88795 56.01787 47.00077 57.90504
## 2017 Q3       66.21274 61.46194 70.96353 58.94702 73.47845
## 2017 Q4       71.23205 65.85721 76.60690 63.01194 79.45217

fc1 %>% checkresiduals()

## 
##  Ljung-Box test
## 
## data:  Residuals from Holt-Winters' additive method
## Q* = 8.6117, df = 3, p-value = 0.03493
## 
## Model df: 8.   Total lags used: 11

fc2 %>% checkresiduals()

## 
##  Ljung-Box test
## 
## data:  Residuals from Holt-Winters' multiplicative method
## Q* = 6.1188, df = 3, p-value = 0.106
## 
## Model df: 8.   Total lags used: 11

fc1$model %>% autoplot()

fc2$model %>% autoplot()

Holt Winters with Dampening

F Gardner uses mostly this

So does Rob if he can only use 1 tool

If you only have 1 tool this is it

Lab Session 8

fc1 <- mytimeseries %>% hw(h=12, seasonal='multiplicative', damped=FALSE)
fc2 <- mytimeseries %>% hw(h=12, seasonal='multiplicative', damped=TRUE)

fc1 %>% autoplot()

fc2 %>% autoplot()

fc1$model %>% autoplot()

fc2$model %>% autoplot()

fc1 %>% summary()

## 
## Forecast method: Holt-Winters' multiplicative method
## 
## Model Information:
## Holt-Winters' multiplicative method 
## 
## Call:
##  hw(y = ., h = 12, seasonal = "multiplicative", damped = FALSE) 
## 
##   Smoothing parameters:
##     alpha = 0.6979 
##     beta  = 0.0087 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 63.7459 
##     b = 1.0282 
##     s=0.9921 0.933 0.9936 1.2091 1.013 1.0172
##            0.979 0.995 0.9791 0.9381 0.9746 0.9761
## 
##   sigma:  0.0468
## 
##      AIC     AICc      BIC 
## 4414.235 4415.713 4483.398 
## 
## Error measures:
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -0.2054187 9.280748 6.244629 -0.2790438 3.338901 0.3336358
##                    ACF1
## Training set 0.05182697
## 
## Forecasts:
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Apr 2018       245.7755 231.0502 260.5008 223.2551 268.2959
## May 2018       245.7272 227.7016 263.7529 218.1593 273.2952
## Jun 2018       236.8256 216.7102 256.9411 206.0617 267.5896
## Jul 2018       247.4718 223.8776 271.0659 211.3876 283.5559
## Aug 2018       251.7887 225.3713 278.2060 211.3868 292.1905
## Sep 2018       248.0573 219.8081 276.3064 204.8539 291.2606
## Oct 2018       258.0681 226.4893 289.6470 209.7724 306.3638
## Nov 2018       257.3264 223.7535 290.8993 205.9811 308.6717
## Dec 2018       307.4929 264.9803 350.0056 242.4754 372.5105
## Jan 2019       253.0461 216.1567 289.9354 196.6287 309.4634
## Feb 2019       237.8969 201.4787 274.3151 182.2000 293.5938
## Mar 2019       253.2718 212.6985 293.8451 191.2203 315.3233

fc2 %>% summary()

## 
## Forecast method: Damped Holt-Winters' multiplicative method
## 
## Model Information:
## Damped Holt-Winters' multiplicative method 
## 
## Call:
##  hw(y = ., h = 12, seasonal = "multiplicative", damped = TRUE) 
## 
##   Smoothing parameters:
##     alpha = 0.7888 
##     beta  = 1e-04 
##     gamma = 1e-04 
##     phi   = 0.9738 
## 
##   Initial states:
##     l = 63.9173 
##     b = 0.2053 
##     s=0.9948 0.935 0.9961 1.2126 1.0159 1.0202
##            0.9781 0.992 0.9751 0.9339 0.9713 0.975
## 
##   sigma:  0.0472
## 
##      AIC     AICc      BIC 
## 4418.242 4419.898 4491.474 
## 
## Error measures:
##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.5050393 9.194427 6.128123 0.2458937 3.290157 0.3274111
##                     ACF1
## Training set -0.04640084
## 
## Forecasts:
##          Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
## Apr 2018       245.8588 230.9987 260.7189 223.1323 268.5854
## May 2018       244.9277 226.0636 263.7917 216.0775 273.7778
## Jun 2018       235.5175 214.1719 256.8632 202.8722 268.1629
## Jul 2018       245.8983 220.7007 271.0959 207.3620 284.4347
## Aug 2018       250.1604 221.8689 278.4519 206.8922 293.4286
## Sep 2018       246.6482 216.3561 276.9404 200.3204 292.9760
## Oct 2018       257.2616 223.3454 291.1778 205.3912 309.1320
## Nov 2018       256.1870 220.2465 292.1275 201.2208 311.1532
## Dec 2018       305.7881 260.4472 351.1290 236.4452 375.1310
## Jan 2019       251.2252 212.0682 290.3821 191.3398 311.1106
## Feb 2019       235.7881 197.3285 274.2477 176.9692 294.6070
## Mar 2019       250.8599 208.1987 293.5212 185.6152 316.1047

fc1 %>% checkresiduals()

## 
##  Ljung-Box test
## 
## data:  Residuals from Holt-Winters' multiplicative method
## Q* = 41.483, df = 8, p-value = 1.693e-06
## 
## Model df: 16.   Total lags used: 24

fc2 %>% checkresiduals()

## 
##  Ljung-Box test
## 
## data:  Residuals from Damped Holt-Winters' multiplicative method
## Q* = 42.336, df = 7, p-value = 4.479e-07
## 
## Model df: 17.   Total lags used: 24

fc1 %>% accuracy()

##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set -0.2054187 9.280748 6.244629 -0.2790438 3.338901 0.3336358
##                    ACF1
## Training set 0.05182697

fc2 %>% accuracy()

##                     ME     RMSE      MAE       MPE     MAPE      MASE
## Training set 0.5050393 9.194427 6.128123 0.2458937 3.290157 0.3274111
##                     ACF1
## Training set -0.04640084

ETS (exponential smoothing)

Exponential smoothing methods: * return point forecasts

Innovations state space models: * general pobabalistic stuff

estimating ETS models:

ets(mytimeseries)

## ETS(M,Ad,M) 
## 
## Call:
##  ets(y = mytimeseries) 
## 
##   Smoothing parameters:
##     alpha = 0.7718 
##     beta  = 0.0025 
##     gamma = 0.0561 
##     phi   = 0.98 
## 
##   Initial states:
##     l = 64.0899 
##     b = 0.2813 
##     s=0.9869 0.9362 1.0043 1.1705 1.0216 1.0293
##            0.9848 0.9998 0.993 0.9436 0.9694 0.9606
## 
##   sigma:  0.0464
## 
##      AIC     AICc      BIC 
## 4404.041 4405.697 4477.272

ets(mytimeseries) %>% autoplot()

ets(mytimeseries) %>% forecast() %>% autoplot()

USAccDeaths %>% 
  ets %>% 
  forecast ->
forecasted_stuff 

head(forecasted_stuff)

## $model
## ETS(A,N,A) 
## 
## Call:
##  ets(y = .) 
## 
##   Smoothing parameters:
##     alpha = 0.5946 
##     gamma = 0.002 
## 
##   Initial states:
##     l = 9248.3628 
##     s=-51.3449 -255.3528 218.2901 -121.771 970.7387 1683.237
##            756.092 306.4212 -489.5627 -739.9004 -1537.792 -739.0552
## 
##   sigma:  292.6907
## 
##      AIC     AICc      BIC 
## 1140.145 1148.716 1174.295 
## 
## $mean
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 1979  8397.497  7599.221  8396.595  8646.510  9443.412  9893.065 10819.124
## 1980  8397.497  7599.221  8396.595  8646.510  9443.412  9893.065 10819.124
##            Aug       Sep       Oct       Nov       Dec
## 1979 10107.105  9015.233  9354.407  8880.418  9085.492
## 1980 10107.105  9015.233  9354.407  8880.418  9085.492
## 
## $level
## [1] 80 95
## 
## $x
##        Jan   Feb   Mar   Apr   May   Jun   Jul   Aug   Sep   Oct   Nov
## 1973  9007  8106  8928  9137 10017 10826 11317 10744  9713  9938  9161
## 1974  7750  6981  8038  8422  8714  9512 10120  9823  8743  9129  8710
## 1975  8162  7306  8124  7870  9387  9556 10093  9620  8285  8466  8160
## 1976  7717  7461  7767  7925  8623  8945 10078  9179  8037  8488  7874
## 1977  7792  6957  7726  8106  8890  9299 10625  9302  8314  8850  8265
## 1978  7836  6892  7791  8192  9115  9434 10484  9827  9110  9070  8633
##        Dec
## 1973  8927
## 1974  8680
## 1975  8034
## 1976  8647
## 1977  8796
## 1978  9240
## 
## $upper
##                80%       95%
## Jan 1979  8772.595  8971.160
## Feb 1979  8035.616  8266.630
## Mar 1979  8886.679  9146.115
## Apr 1979  9184.957  9469.994
## May 1979 10026.222 10334.743
## Jun 1979 10517.092 10847.432
## Jul 1979 11481.809 11832.614
## Aug 1979 10806.314 11176.454
## Sep 1979  9749.151 10137.664
## Oct 1979 10121.466 10527.522
## Nov 1979  9679.242 10102.115
## Dec 1979  9915.072 10354.225
## Jan 1980  9256.534  9711.281
## Feb 1980  8486.738  8956.562
## Mar 1980  9311.707  9796.138
## Apr 1980  9588.408 10087.018
## May 1980 10411.355 10923.753
## Jun 1980 10886.371 11412.195
## Jul 1980 11837.160 12376.076
## Aug 1980 11149.285 11700.983
## Sep 1980 10081.011 10645.200
## Oct 1980 10443.272 11019.681
## Nov 1980  9991.889 10580.266
## Dec 1980 10219.269 10819.454
## 
## $lower
##                80%      95%
## Jan 1979  8022.399 7823.834
## Feb 1979  7162.825 6931.812
## Mar 1979  7906.510 7647.075
## Apr 1979  8108.063 7823.026
## May 1979  8860.602 8552.081
## Jun 1979  9269.038 8938.698
## Jul 1979 10156.438 9805.634
## Aug 1979  9407.895 9037.756
## Sep 1979  8281.314 7892.801
## Oct 1979  8587.349 8181.293
## Nov 1979  8081.593 7658.721
## Dec 1979  8255.912 7816.759
## Jan 1980  7538.459 7083.713
## Feb 1980  6711.703 6241.880
## Mar 1980  7481.483 6997.052
## Apr 1980  7704.612 7206.001
## May 1980  8475.469 7963.070
## Jun 1980  8899.759 8373.935
## Jul 1980  9801.087 9262.171
## Aug 1980  9064.924 8513.227
## Sep 1980  7949.455 7385.266
## Oct 1980  8265.543 7689.133
## Nov 1980  7768.947 7180.570
## Dec 1980  7951.715 7351.530

forecasted_stuff %>%
  autoplot

h02 %>% ets %>% forecast %>% autoplot

h02 %>% ets %>% print

## ETS(M,Ad,M) 
## 
## Call:
##  ets(y = .) 
## 
##   Smoothing parameters:
##     alpha = 0.1953 
##     beta  = 1e-04 
##     gamma = 1e-04 
##     phi   = 0.9798 
## 
##   Initial states:
##     l = 0.3945 
##     b = 0.0085 
##     s=0.874 0.8197 0.7644 0.7693 0.6941 1.2838
##            1.326 1.1765 1.1621 1.0955 1.0422 0.9924
## 
##   sigma:  0.0676
## 
##        AIC       AICc        BIC 
## -122.90601 -119.20871  -63.17985

autoplot(a10)

a10 %>% ets %>% autoplot

jdl_model <- function(dataset) {
  dataset %>% ets -> fit_ets
  print(fit_ets)
  print(fit_ets %>% autoplot)
  print(fit_ets %>% forecast %>% autoplot)
}

jdl_model(bicoal)

## ETS(M,N,N) 
## 
## Call:
##  ets(y = .) 
## 
##   Smoothing parameters:
##     alpha = 0.8205 
## 
##   Initial states:
##     l = 542.665 
## 
##   sigma:  0.1262
## 
##      AIC     AICc      BIC 
## 595.2499 595.7832 600.9253

jdl_model(chicken)

## ETS(M,N,N) 
## 
## Call:
##  ets(y = .) 
## 
##   Smoothing parameters:
##     alpha = 0.98 
## 
##   Initial states:
##     l = 159.8322 
## 
##   sigma:  0.1691
## 
##      AIC     AICc      BIC 
## 635.2382 635.6018 641.9836

jdl_model(dole)

## ETS(M,Ad,M) 
## 
## Call:
##  ets(y = .) 
## 
##   Smoothing parameters:
##     alpha = 0.697 
##     beta  = 0.124 
##     gamma = 0.303 
##     phi   = 0.902 
## 
##   Initial states:
##     l = 2708.6621 
##     b = 836.017 
##     s=1.0404 0.8893 0.9103 1.0301 1.0576 1.0584
##            0.9801 0.9632 1.021 0.9838 1.0145 1.0514
## 
##   sigma:  0.0935
## 
##      AIC     AICc      BIC 
## 10602.67 10604.30 10676.19

jdl_model(usdeaths)

## ETS(A,N,A) 
## 
## Call:
##  ets(y = .) 
## 
##   Smoothing parameters:
##     alpha = 0.5972 
##     gamma = 0.0019 
## 
##   Initial states:
##     l = 9195.6403 
##     s=-62.6129 -270.0351 263.3823 -89.4907 1005.529 1662.647
##            795.2585 333.326 -551.161 -737.5102 -1552.872 -796.4611
## 
##   sigma:  294.4663
## 
##      AIC     AICc      BIC 
## 1141.016 1149.587 1175.166

jdl_model(bricksq)

## ETS(M,A,M) 
## 
## Call:
##  ets(y = .) 
## 
##   Smoothing parameters:
##     alpha = 0.973 
##     beta  = 0.0032 
##     gamma = 1e-04 
## 
##   Initial states:
##     l = 196.0437 
##     b = 4.9839 
##     s=1.0041 1.0598 1.0269 0.9093
## 
##   sigma:  0.0514
## 
##      AIC     AICc      BIC 
## 1725.473 1726.714 1752.864

jdl_model(lynx)

## ETS(M,N,N) 
## 
## Call:
##  ets(y = .) 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
## 
##   Initial states:
##     l = 2372.8047 
## 
##   sigma:  0.9594
## 
##      AIC     AICc      BIC 
## 2058.138 2058.356 2066.346

jdl_model(ibmclose)

## ETS(A,N,N) 
## 
## Call:
##  ets(y = .) 
## 
##   Smoothing parameters:
##     alpha = 0.9999 
## 
##   Initial states:
##     l = 459.9339 
## 
##   sigma:  7.2637
## 
##      AIC     AICc      BIC 
## 3648.450 3648.515 3660.182

jdl_model(eggs)

## ETS(M,N,N) 
## 
## Call:
##  ets(y = .) 
## 
##   Smoothing parameters:
##     alpha = 0.8198 
## 
##   Initial states:
##     l = 278.8889 
## 
##   sigma:  0.1355
## 
##      AIC     AICc      BIC 
## 1043.286 1043.553 1050.916

jdl_model(ausbeer)

## ETS(M,A,M) 
## 
## Call:
##  ets(y = .) 
## 
##   Smoothing parameters:
##     alpha = 0.2087 
##     beta  = 0.0302 
##     gamma = 0.1984 
## 
##   Initial states:
##     l = 255.892 
##     b = 0.6795 
##     s=1.183 0.9113 0.8614 1.0444
## 
##   sigma:  0.0361
## 
##      AIC     AICc      BIC 
## 2354.249 2355.114 2384.709

## works but needs to print name on every run of the function

data_list <- list(bicoal, chicken, dole, usdeaths, bricksq, lynx, ibmclose, eggs, ausbeer)
data_list$names <- c("bicoal", "chicken", "dole", "usdeaths", "bricksq", "lynx", "ibmclose", "eggs", "ausbeer")
lapply(data_list, jdl_model)

autoplot(mytimeseries)

train <- window(mytimeseries, end=c(2010,12))
test  <- window(mytimeseries, start=2011)

f1 <- snaive(train, h=length(test))
f2 <- rwf(train, h=length(test))
f3 <- rwf(train, drift = TRUE, h=length(test))
f4 <- meanf(train, h=length(test))
f5 <- hw(train, h=length(test), seasonal = 'multiplicative')
f6 <- ets(train) %>% forecast(h=length(test))

c(
  SN=accuracy(f1, test)[2,"RMSE"],
  rwf1=accuracy(f2, test)[2,"RMSE"],
  rwf2=accuracy(f3, test)[2,"RMSE"],
  meanf=accuracy(f4, test)[2,"RMSE"],
  hw=accuracy(f5, test)[2,"RMSE"],
  ets=accuracy(f6, test)[2,"RMSE"]
)

##       SN     rwf1     rwf2    meanf       hw      ets 
## 23.40458 37.09060 56.64803 58.73867 59.85277 25.43061

now let’s do time series cross validation

e1 <- tsCV(mytimeseries, snaive, h=12)
e2 <- tsCV(mytimeseries, naive, h=12)
e3 <- tsCV(mytimeseries, rwf, drift=TRUE, h=12)
e4 <- tsCV(mytimeseries, meanf, h=12)
e5 <- tsCV(mytimeseries, hw, h=12, seasonal='multiplicative')

etsfc <- function(y,h){
  y %>% ets(model="MAM", damped=TRUE) %>%
    forecast(h=h)
}

e6 <- tsCV(mytimeseries, etsfc, h=12)

MSE <- cbind(
  h=1:12,
  SNaive=colMeans(tail(e1, -14)^2, na.rm=TRUE), 
  Naive=colMeans(tail(e2, -14)^2, na.rm=TRUE), 
  drift=colMeans(tail(e3, -14)^2, na.rm=TRUE), 
  Mean=colMeans(tail(e4, -14)^2, na.rm=TRUE), 
  HW=colMeans(tail(e5, -14)^2, na.rm=TRUE), 
  ETS=colMeans(tail(e6, -14)^2, na.rm=TRUE)
)
MSE

##       h   SNaive     Naive     drift     Mean       HW      ETS
## h=1   1 695.0202  411.9992  413.3894 4708.335 106.7760 104.3003
## h=2   2 696.6544  543.7812  547.3381 4757.235 171.8868 158.4553
## h=3   3 698.2653  596.3238  601.7417 4806.157 251.5070 225.1077
## h=4   4 699.9288  661.9810  668.9494 4855.374 323.9946 294.1775
## h=5   5 701.6235  755.7454  765.1886 4904.751 398.8013 364.3042
## h=6   6 703.3177  917.1529  931.6504 4954.304 479.8465 418.0667
## h=7   7 705.0283  878.4346  893.1175 5003.275 559.1300 494.1943
## h=8   8 706.7069  905.6455  922.2475 5053.316 653.4687 555.8380
## h=9   9 708.3497  950.4590  969.6496 5103.750 730.6287 624.1344
## h=10 10 710.0246 1010.2253 1032.6967 5154.446 822.1079 694.4674
## h=11 11 711.6948  983.8066 1005.4918 5205.465 916.4102 757.7232
## h=12 12 713.3165  713.3165  722.5366 5257.032 996.6152 808.3827

day 2

Transformations. Box Cox

elec %>% autoplot

(lambda <- BoxCox.lambda(elec))

## [1] 0.2654076

elec %>%
  BoxCox(lambda) %>%
  autoplot

fc <- rwf(eggs, drift=TRUE, lambda=0, h=50, level=80)
fc2 <- rwf(eggs, drift=TRUE, lambda=0, h=50, level=80,
  biasadj=TRUE)
autoplot(eggs) +
  autolayer(fc, series="Simple back transformation") +
  autolayer(fc2, series="Bias adjusted", PI=FALSE) +
  guides(colour=guide_legend(title="Forecast"))

lambda_adjusted <- function(data) {
  data %>% autoplot

  lambda <- BoxCox.lambda(data)
  print(lambda)
  
  print("with transform")
  data %>%
    BoxCox(lambda) %>%
    forecast %>%
    autoplot
  
  print("without transform")
  data %>%
    forecast %>%
    autoplot
}

data_list <- list(usnetelec,
                  mcopper,
                  enplanements,
                  a10
                  )
lapply(data_list, lambda_adjusted)

## [1] 0.5167714
## [1] "with transform"
## [1] "without transform"
## [1] 0.1919047
## [1] "with transform"
## [1] "without transform"
## [1] -0.2269461
## [1] "with transform"
## [1] "without transform"
## [1] 0.1313326
## [1] "with transform"
## [1] "without transform"

## [[1]]

## 
## [[2]]

## 
## [[3]]

## 
## [[4]]

Lab Sessions: Forecasting using R

JD Long

25 June 2018